Subsurface Oxygen \(O_2\)#

Hide code cell source
import warnings
warnings.filterwarnings("ignore")
import os
import os.path as op
import sys

import pandas as pd
import numpy as np
import xarray as xr
import geopandas as gpd
import cartopy.crs as ccrs
import matplotlib.pyplot as plt

sys.path.append("../../../../indicators_setup")
from ind_setup.plotting_int import plot_timeseries_interactive, plot_oni_index_th
from ind_setup.plotting import plot_base_map, plot_map_subplots, add_oni_cat, plot_bar_probs, fontsize

sys.path.append("../../../functions")
from data_downloaders import download_oni_index

Setup#

Define area of interest

#Area of interest
lon_range  = [129.4088, 137.0541]
lat_range = [1.5214, 11.6587]

EEZ shapefile

path_figs = "../../../matrix_cc/figures"
shp_f = op.join(os.getcwd(), '..', '..','..', 'data/Palau_EEZ/pw_eez_pol_april2022.shp')
shp_eez = gpd.read_file(shp_f)

Load Data#

Superficial depth

data_xr = xr.open_dataset(op.join(os.getcwd(), '..', '..','..', 'data/data_phyc_o2_ph.nc')).isel(depth = 0)
dataset_id = 'o2'
label = 'O2 [µmol/kg]'

Analysis#

Plotting#

Average#

ax = plot_base_map(shp_eez = shp_eez, figsize = [10, 6])
im = ax.pcolor(data_xr.longitude, data_xr.latitude, data_xr.mean(dim='time')[dataset_id], transform=ccrs.PlateCarree(), 
                cmap = 'BuPu', 
                vmin = np.nanpercentile(data_xr.mean(dim = 'time')[dataset_id], 1), 
                vmax = np.nanpercentile(data_xr.mean(dim = 'time')[dataset_id], 99))
ax.set_extent([lon_range[0], lon_range[1], lat_range[0], lat_range[1]], crs=ccrs.PlateCarree())
plt.colorbar(im, ax=ax, label= label)
plt.savefig(op.join(path_figs, 'F17_O2_mean_map.png'), dpi=300, bbox_inches='tight')
../../../_images/9fbfc0da983ed622786bc1aa3cc3f58ea4e4a5b3159522ca6e6bca4fe362df1f.png

Annual average#

data_y = data_xr.resample(time='1YE').mean()
im = plot_map_subplots(data_y, dataset_id, shp_eez = shp_eez, cmap = 'BuPu', 
                  vmin = np.nanpercentile(data_xr.min(dim = 'time')[dataset_id], 1), 
                  vmax = np.nanpercentile(data_xr.max(dim = 'time')[dataset_id], 99),
                  cbar = 1)
../../../_images/7bae6337383efdb7fd289b17c4f57590b4ce599a68354d616f5ce11a8323a00d.png

Annual anomaly#

data_an = data_y - data_xr.mean(dim='time')
fig = plot_map_subplots(data_an, dataset_id, shp_eez = shp_eez, cmap='RdBu_r', vmin=-4, vmax=4, cbar = 1)
../../../_images/2df149e6c58e5057141dea2c59679b5788d9538ab4163e84b70b988a1f7e1fa4.png

Average over area#

dict_plot = [{'data' : data_xr.mean(dim = ['longitude', 'latitude']).to_dataframe(), 
              'var' : dataset_id, 'ax' : 1, 'label' : f'{label} - MEAN AREA'},]
fig = plot_timeseries_interactive(dict_plot, trendline=True, scatter_dict = None, figsize = (25, 12), label_yaxes = label);
fig.write_html(op.join(path_figs, 'F17_O2_mean_trend.html'), include_plotlyjs="cdn")

Top 10 years#

data_mean = data_xr.mean(dim = ['longitude', 'latitude']).to_dataframe()
data_mean = data_mean.resample('YE').mean()
top_10 = data_mean.sort_values(by='o2', ascending=False).head(10)

fig, ax, trend = plot_bar_probs(x=data_mean.index.year, y=data_mean.o2, trendline=True,
                                y_label='Dissolved O2', figsize=[15, 4], return_trend=True)
ax.set_ylim(data_mean.o2.min()-0.5, data_mean.o2.max()+0.5)

im = ax.scatter(top_10.index.year, top_10.o2,
                c=top_10.o2.values, s=100, ec='pink', cmap='rainbow', label='Top 10 warmest years')
plt.colorbar(im, pad=0, shrink=.7).set_label('Absolute Dissolved O2', fontsize=fontsize)
ax.set_title('Annual Mean Dissolved O2', fontsize=15);
../../../_images/d97bb2571faa516a9ef619dc394391fd3e1fefe55a4dcf85e619a45665312653.png

Timeseries at a given point#

loc = [7.37, 134.7]
dict_plot = [{'data' : data_xr.sel(longitude=loc[1], latitude=loc[0], method='nearest').to_dataframe(), 
              'var' : dataset_id, 'ax' : 1, 'label' : f'{label} at [{loc[0]}, {loc[1]}]'},]
ax = plot_base_map(shp_eez = shp_eez, figsize = [10, 6])
ax.set_extent([lon_range[0], lon_range[1], lat_range[0], lat_range[1]], crs=ccrs.PlateCarree())
ax.plot(loc[1], loc[0], '*', markersize = 12, color = 'royalblue', transform=ccrs.PlateCarree(), label = 'Location Analysis')
ax.legend()
<matplotlib.legend.Legend at 0x1828eecc0>
../../../_images/038a2757803511e8b284fa11a57d0ba084269dbabe56e09215299e69572e231d.png
fig = plot_timeseries_interactive(dict_plot, trendline=True, scatter_dict = None, figsize = (25, 12), label_yaxes = label);

ONI index analysis#

p_data = 'https://psl.noaa.gov/data/correlation/oni.data'
df1 = download_oni_index(p_data)
lims = [-.5, .5]
plot_oni_index_th(df1, lims = lims)

Group by ONI category

df1 = add_oni_cat(df1, lims = lims)
df1['ONI'] = df1['oni_cat']
data_xr['ONI'] = (('time'), df1.iloc[np.intersect1d(data_xr.time, df1.index, return_indices=True)[2]].ONI.values)
data_xr['ONI_cat'] = (('time'), np.where(data_xr.ONI < lims[0], -1, np.where(data_xr.ONI > lims[1], 1, 0)))
data_oni = data_xr.groupby('ONI_cat').mean()

Average#

fig = plot_map_subplots(data_oni, dataset_id, shp_eez = shp_eez, cmap = 'BuPu', 
                  vmin = np.nanpercentile(data_xr.mean(dim = 'time')[dataset_id], 1) - 1, 
                  vmax = np.nanpercentile(data_xr.mean(dim = 'time')[dataset_id], 99) + 1,
                  sub_plot= [1, 3], figsize = (20, 9),  cbar = True, cbar_pad = 0.1,
                  titles = ['La Niña', 'Neutral', 'El Niño'],)
plt.savefig(op.join(path_figs, 'F17_O2_ENSO.png'), dpi=300, bbox_inches='tight')
../../../_images/e9e6fd4dc8566874c783df73e47e43059515af14480ce2323c646b166b8a7f74.png

Anomaly#

data_an = data_oni - data_xr.mean(dim='time')
fig = plot_map_subplots(data_an, dataset_id, shp_eez = shp_eez, cmap='RdBu_r', vmin=-1, vmax=1,
                  sub_plot= [1, 3], figsize = (20, 9),  cbar = True, cbar_pad = 0.1,
                  titles = ['La Niña', 'Neutral', 'El Niño'],)
../../../_images/e0ee38d69bc62511cc190fbe37b141917b211ff22add3f3ebbd5148b5ee120a8.png